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Abstract 

We present phase-space techniques for the modelling of spontaneous emission in 
two-level bosonic atoms. The positive-P representation is shown to give a full and 
complete description and can be further developed to give exact treatments of the 
interaction of degenerate bosons with the electromagnetic field in a given experi- 
mental situation. The Wigner representation, even when truncated at second order, 
is shown to need a doubling of the phase-space to allow for a positive-definite diffu- 
sion matrix in the appropriate Fokker-Planck equation and still fails to agree with 
the full quantum results of the positive-P representation. We show that quantum 
statistics and correlations between the ground and excited states affect the dynamics 
of the emission process, so that it is in general non-exponential. 
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1 Introduction 



The study of spontaneous emission has a long history, beginning with the fa- 
mous Einstein A and B coefficients [1] and the Wigner- Weisskopf law [2]. A 
comprehensive account can be found in Agarwal [3]. Previous studies have con- 
sidered atoms which are not within one de Broglie wavelength of each other, 
so that the fermionic or bosonic nature of the atoms need not be considered. 
This condition holds for dilute high-temperature atomic samples, but not for 
degenerate quantum gases. In particular, it does not hold for Bose-Einstein 
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condensates [4] and now that these are readily available in laboratories around 
the world, it is of some interest to develop methods which will allow for a full 
quantum treatment of both the atomic and electromagnetic fields. In this pa- 
per we develop and demonstrate phase-space methods, using the positive-P [5] 
and Wigner [6] representations, for the treatment of spontaneous emission from 
two-level bosonic atoms. We use a simple model without spatial dependence 
of the atomic fields and without collisional and dipole-dipole interactions, al- 
though these can readily be included for larger atomic samples, along with 
other details which would be necessary to model a realistic experimental sit- 
uation. The simplicity of our model allows us to consider different quantum 
states of the atomic fields and demonstrate clearly that correlations which can 
build up between the ground and excited states have a noticeable effect on the 
dynamics. We can also give a clear demonstration of the difference between 
the quantum mechanical predictions of the positive-P representation and an 
approximation which we derive from the Wigner representation. 



2 Formalism 

To begin, we define ip(x) as a wavefunction for a bosonic atomic field, with 
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the vfij (x) are then wave-functions for different (distinguishable) atomic fields. 
In second quantisation, the expansion coefficients are changed to operators 
so that OjCij becomes the number operator for atoms of the jth type. For 
simplicity, we consider one chemical species of atom with two electronic levels, 
j = a, b, for the ground and excited states, respectively, and set a a = a and 
a b = b. Having introduced the second-quantised picture we could also include 
fermionic atoms, for which Eq. 1 is not defined, and use fermionic phase- 
space methods which are under development [7] . However, in this work we will 
consider only bosonic atoms. Since the total number of atoms is conserved, 
the physically relevant (bosonic) operator combinations can be defined by 



a ] b\n a , n b ) = \J (n a + l)n h \n a + l,n b -l), 

&a\n a ,n b ) = ^(n b + l)n a \n a -l,n b + 1), (2) 

where \n a ,n b ) signifies an atomic field with n a atoms in the ground state and 
rib atoms in the excited state. 

We will consider first the case of a single-atom Fock state, where (a^a+b>b) = 1. 
The atomic Fock space then reduces to two dimensions, with basis vectors 
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|0, 1) and 1 1,0). The combinations of operators (2) then play an analogous 
role to the normal spin raising and lowering operators. We can now make the 
following correspondences with the normal Pauli spin operators: 

a + <-> Pa, a~ <-> a)b, o z <-> ^ (Pb — a^a) , (3) 
and look at the commutation relations. We find 



Pa, a)b = Pb — a) a, 

= -2Pa, 
= 2a r b, 



Pa, Pb — a) a 
a r b, Pb 



a) a 



(4) 



all completely equivalent to those of two-level atomic operators such as used 
in the well-known Jaynes-Cummings model [8], as long as we are considering 
a single atom. We note here that a similar formalism has previously been used 
by Bonifacio and Preparata [9], although they did not develop a phase-space 
representation of the problem as we will do in what follows. 



Not unexpectedly, this simple relation to the operators used in the Jaynes- 
Cummings model does not survive beyond the single-atom Fock-state case. 
Consider for example the expressions for the probability of an atom being in 
a particular state. In the latter, (with Pjj signifying the probability of being 
in state j) 

P bb =(a+a-}, P aa = (a-a + ), (5) 
and, if we naively use these correspondences (3), we find 



+ <-> (Patfb) = (Pb) + (Pba ] a), 



(a~a + ) <-> (aWa) = (a) a) + (Pba)a), (6) 

which are obviously not the correct probabilities. This is because the Pauli 
operators describe a single atom with one fermion (electron) which can be in 
either the ground or excited state. This means that the density matrix need 
only be 2 x 2, whereas we wish to investigate bosonic atoms where the density 
matrix is, in priciple, infinite. To solve this problem, we define 

p - @h P (?) 

bb (Pb) + (ata) ' " (Pb) + (ata) ' 1 ' 

which are now the correct probabilities. We may define the atomic coherences 
as 

P ba = (a-)^(a^b), P ah = (a + )^(Pa). (8) 
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We may now write any master equation in terms of these bosonic operators 
and use the standard mappings [10] to find Fokker-Planck and stochastic dif- 
ferential equations in the relevant phase-space representations. As an example 
to demonstrate our method, we will examine spontaneous emission from ex- 
cited bosonic atoms. 



3 Spontaneous decay 



Spontaneous decay into a zero temperature thermal bath can be modelled by 
the Hamiltonian 

H bath = h (aPv + tfbT*) , (9) 

where the V are operators for the bath quanta. Following the usual meth- 
ods [11], this leads to the master equation 

^ = | (2a%ati - dPa% - patitfb) . (10) 



Looking at Eq. 10, we can immediately find equations for expectation values 
as 



dt 1 dt 

We find 



d{£) -Tr(M. (U) 



*^ = - K (@b) + (tfcftb)), 

(Jj V 

d(a)d) 



dt 



k (<S + 6> + (tfaftb)) , (12) 



so that any correlations which build up between the two atomic fields can be 
important in the dynamics of the decay. These may be expected to depend on 
the initial atomic quantum states, something which is easily modelled using 
positive-P and Wigner representations. In the standard techniques which use 
a density matrix approach, approximations are usually necessary to model 
anything but number states, which are not the only possible choice when we 
wish to treat degenerate quantum gases. 

The choice of Hamiltonian (9) shows clearly that our analysis will not de- 
scribe superrradiant emission. The latter may be defined as a feedback of the 
radiated field on atoms, whereas, by using a bath for the emitted field this 
feedback is formally eliminated. Irrespective of how important superrradiation 
may happen to be in a real experiment, the effects we describe in this paper 
are rooted in the many-body atom-field statistics. 
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4 Positive-P representation 



We will first develop our stochastic equations in the positive-P representation. 
Using the usual techniques [10], the master equation (10) is mapped onto a 
Fokker-Planck equation for the P-function [12,13], 



dP 
~dt 



K 

2 



d 



d 



d 



^\Pfa + ^-\(3\ 2 a*-i-(\a\ 2 + l)(3 



d 



+- 



da 
d 2 



da 



dp 



(\a\ 2 + l)(3* 



dada 



-m 2 - 



d 2 



dad/3 



af3 



d 2 



da*d(3* 



* ft* 

a p 



d[3 

P(a,P,t). (13) 



The diffusion matrix, D, of Eq. 13, is not positive-definite so we must use the 
positive-P representation. The matrix can be factorised as 



B = \[k 



where D 
equations 



Soft 





;Va^ 



\\fa~P 



\\fa~P 







^^+P+ ±Va~+p 



+ 













(14) 



BB T . This immediately allows us to write the Ito stochastic 



K n+n 1 I n , ■ \ / ^P^P / . n 

-P^pa - -^naP (r} 3 - vq{) + W — ^— {rj 5 + ir] & ) , 
^P + pa + - ^Ka+p+ (774 - irj 2 ) + ]J ~*£Jl ( m _ i m ) , 
~ (a + a + l) p + ^\JnaP (r) 3 + ir^) , 

~ (a + a + l)p + + y*a+P+ (r? 4 + iffr) , (15) 

which may be numerically integrated. Note that, as always in the positive-P 
representation, there is a correspondence between the operators a, b, a\ W and 
the c-number variables a, (3, a + , P + , such that 

{a+) m {P+) n avp<i -> ((a t ) m (S t )'W), (16) 

where the left hand side is a classical average over stochastic trajectories and 
the right hand side is a quantum-mechanical expectation value. 

Before we resort to stochastic integration, we will see what information we can 
get from the Ito equations about the average decay rates, with the caveat that 



da 
~dt 
da + 
~dT 

dp 

~dt 

dp+ 

~df 
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a mean-field approach may not be particularly meaningful for small atomic 
samples. Using the rules of Ito calculus, where 

d(xy) = y dx + x dy + dy dx, (17) 

we find mean-field equations for the populations in the ground (N a = a + a) 
and excited (N b = f3 + f3) states, 



^ = n(N a + l)N b , 
dN b 



dt 



= -K(N a + l)N b . (18) 



For those not familiar with stochastic calculus, we note that merely dropping 
the noise terms from Eq. 15 would not give the correct equations for the 
populations, although this would be the correct procedure if we were using 
Stratonovich equations. With the total number of atoms, we may now 
write an equation for N b , 



^± = - K (N a + l)N b 

— —k (Nt — N b + 1) N b , (19) 

which may be solved to give 

at (t) = ?m i^i ± *) (20) 

b{) N b (0) [N T + 1- N b (0)} e (^+i)«t • [ZV) 

In the case of one atom initially excited, this simplifies to 

N b (t) = sech(Kt)e~ K \ (21) 

although we would not expect mean-field solutions to be accurate for such 
small numbers of quanta. However, we do see that the decay rate is propor- 
tional to the number of quanta which can be in the final state, in a manner 
reminiscent of Fermi's golden rule [14]. A work on superradiance, by Rehler 
and Eberly [15], gives an expression for the atomic energy (their Eq. 5.1), 
which, in the limit that all atoms are within a wavelength of each other so 
that their parameter jj, is equal to 1, becomes 

W{t) = ~N b {0) [e^°) +1 ) Kt - (N b (0) + 2)] / ( e ^^ +1 ^* + N b (0)) . (22) 

We note that, once we redefine their energy scale so that the energy of the 
ground state is zero rather than negative and all atoms are initially excited, 
this gives the same result as our Eq. 20. (In Ref. [15], the energy of a ground 
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(excited) state atom is defined as —(+)hu/2.) It is interesting that the same 
result has been arrived at by different means, although our result allows for 
some initial population in the ground state and, as shown below by the numer- 
ical results, is only an approximation to the real result, which we will show to 
depend on the quantum state of the atomic ensemble. In the regime where all 
atoms were within an optical wavelength, we would also expect dipole-dipole 
forces and maybe even collisional processes to become important. We stress 
here that the non-exponential decay predicted by our approach is not the clas- 
sical superradiance, but rather is due to correlations between the ground and 
excited state atomic fields. 



To obtain the full quantum results, the Stratonovich version of Eq. 15 was nu- 
merically integrated. As usual in quantum stochastic integration, we must 
decide on which initial conditions to use for our atomic fields. The sim- 
plest initial condition to model in the positive-P representation is a coherent 
state, \cxo), which can be represented by the pseudoprobability distribution 
P(a, a + ) = 5{a — ao)5(a + — ccq). Other quantum states which may arise nat- 
urally when we treat atoms are the chaotic state and the Fock state, which 
has a fixed number of atoms. It has been shown that any quantum state may 
be represented by the positive-P distribution [5] 



P(a,a + ) 



4tt 2 



\ 2 



a + (a + Y 



(23) 



which we wish to sample for a Fock state with p = \n)(n\. Introducing the 
new variables 



a — (a 



H = 

we find the separable expression 



and 7 



a + («+)* 



(24) 



-M 2 U|2« e -M 2 



TT 7TO! 

e -w 2 r(| 7 |V 



71 71 

where 



(25) 



r( *' n) = (^ (26) 

is the Gamma distribution. The variable /x is easily sampled via standard 
methods, while the Gamma distribution is sampled using a method given by 
Marsaglia and Tsang [16] to give z = |7| 2 , so that 7 = ^fz e 10 , where 6 is 
uniform on [0, 2ir). We then invert to find 

a = ji + 7 and a + = 7* — //, (27) 

which are now correctly distributed to represent the Fock state \n). For n — 1, 
the decay follows the well-known exponential law, as asking whether a single 
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Fig. 1. Positive-P solution averaged over 3.4 x 10 5 stochastic trajectories (solid 
line) for decay of an N = 1 Fock state, with k = 0.4. The exponential solution 
(dash-dotted line) and the mean-field analytical solution (dotted line) are shown for 
comparison. 

isolated atom is bosonic or fermionic is a meaningless question. In fact, it is 
much simpler to calculate the decay for Fock states using a master equation 
formalism rather than stochastic equations [17] and we present the result 
shown in Fig. 1 more as evidence that our phase-space approach gives reliable 
results. We find that the stochastic results for larger n agree with the master 
equation results wherever the integration converges. This gives us confidence 
that the results for other quantum states, for which a master equation solution 
is not so easily found, will also be accurate. 

In Fig. 2 we show results for an initial coherent state, |/3), where (3 = 1 and, 
in Fig. 3, the correlation which builds up between the ground and excited 
state populations for this initial state. For the initial n = 1 Fock state, this 
correlation was zero. It can be seen that the decays are noticeably different, 
although we again note that as soon as we have the possibility of more than 
one atom being present, collisions and dipole-dipole interactions would also 
play a role. These have not been included in our analysis at this stage as we 
are more interested in a proof of principle rather than modelling exactly a 
given physical system. 

An interesting theoretical application of this model can be made for bosons 
in lattice wells. The Heisenberg uncertainty principle has been used to infer 
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Fig. 2. Stochastic result (solid line) from Eq. 15, with k = 0.2 and A^(0) = 1, 
with the atom initially excited. For comparison, the standard exponential decay 
(dash-dotted line) and the solution from Eq. 20 (dotted line) are also shown. Note 
that in this result and that shown in Fig. 3, the excited atom is initially in a coherent 
state. The equations were integrated over 9.9 x 10 5 trajectories. 

that atoms condensed at lattice sites were in number states, due to an increase 
in phase noise [18]. We note here that number states are not the only states 
with phase noise above the coherent state level, as chaotic and thermal states, 
among others, will also exhibit this property. As small numbers of atoms in a 
lattice site may be a perfect candidate for measurements of the spontaneous 
emission rate, we will investigate this process for a chaotic state. These states 
have a particularly simple P-function, with 



where n is the average number present in the mode [11]. If the state is a 
mixture of coherent and chaotic states, i.e. a chaotic state with a coherent 
displacement, the P-function is written as 



where f3 is the coherent displacement. We show the results of an averaging of 
1.56 x 10 5 stochastic trajectories in Fig. 4, for /3 = 1 and n = 0.1. In Fig. 5, 
we show the results for /3 = and n = 1. By comparison with Fig. 2, we see 
that the decay in the chaotic case is noticeably faster. 






(29) 
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Fig. 3. Stochastic solution for (a'ab'b), for the same parameters as used in Fig. 2. 




Kt 

Fig. 4. Decay of an initially excited ensemble with 1 atom in the coherent component 
and an average of 0.1 in the chaotic component. The stochastic average is shown 
as the solid line, while the analytical solution of Eq. 20 is shown for comparison 
(dotted line), along with an exponential solutions with decay rate k (dash-dotted 
line). 
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Fig. 5. Decay of an initially excited ensemble with n = 1 and (3q = 0. The stochastic 
average is shown as the solid line, while the analytical solution of Eq. 20 is shown 
for comparison (dotted line), along with an exponential solutions with decay rate k 
(dash-dotted line). 



5 Wigner representation 



In the theory of stochastic electrodynamics [19], spontaneous emission is often 
claimed as being explicable as stimulated emission which is actually stimulated 
by vacuum fluctuations [20]. The Wigner representation commonly used in 
quantum optics is equivalent to this theory if the Fokker-Planck equation for 
the Wigner function is truncated at second order, a positive-definite diffusion 
matrix is found, and distributions with positive Wigner functions are con- 
sidered. Using the standard correspondences for Ito calculus [10], the master 
equation (Eq. 10) can be mapped onto a generalised Fokker-Planck equation 
for the Wigner function, 
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dt 
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{- 
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+ 




dp 



d_ 

da 



(dfl 2 - !>«) + £ (d^ 2 - 5K) 



+ - 



9 2 



dada 
d 2 , 



dadp 



2a/?) 



d 2 



d 2 



da* 8(3 



d 3 



d(3dj3 
: (-2a* j?) 
d 3 



(2\a\ 2 + l) 



dada*df3 V 2 / dada* dp* 
dadpdp* \2 J da*dpdp* \2 



W(a,P,t). 



(30) 



Although methods exist for developing stochastic difference equations for a 
system which gives derivatives of higher than second-order in the Fokker- 
Planck equation [21,22], we will truncate the above equation at second order as 
our aim is to compare the predictions of the positive-P representation with that 
of the representation which results from discarding the third-order derivatives 
in Eq. 30. However, what we do notice is that the diffusion matrix of Eq. 30 is 
not positive-definite and thus has no straightforward mapping onto stochastic 
differential equations. As shown in Ref. [22], we may double the phase-space 
in a way analogous to that used in the positive-P representation and find 
stochastic equations for four independent variables in what we may call a 
truncated positive- Wigner representation. 

We therefore make the changes a* — * a + and /3* — > P + , noting that averages of 
these are equivalent to symmetrically ordered operator expectation values, so 
that these are not the same as the variables used in the positive-P equations, 
although this should be obvious by context. One possible factorisation of the 
diffusion matrix of Eq. 30 is then 



B 



w 



A 4 + 4 C 



(31) 



where 



A 



y K (p+p-\) f^ K (p+p-i) 

lyfiT( 










K[a^a 



+ 1) 








+ |) f^a+a + i) 



n{a + a + 



(32) 
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and 



\\J KQ.fi 












C 







iVna+fi+ 



yna+fi+ 



(33) 



\\J KO.fi 












^Ka+fi + 







yKa+fi+ 



and 4 is the 4x4 null matrix. This allows us to write the following stochastic 
differential equations, 



These can again be readily integrated numerically, at least for states which 
have a well-behaved Wigner function. 

We present a numerical solution to this equation, for an initial coherent state 
with an average of one excited atom, in Fig. 6, along with the positive-P solu- 
tion and the mean-field solution of Eq. 20. What we see is that, even though 
the phase-space has been doubled so that we are representing nonclassical 
dynamics, the truncated positive Wigner solution is closer to the mean-field 
prediction than it is to the full quantum solution. If the naive procedure of 
dropping the noise terms completely from Eq. 34 and integrating the resulting 
equations is followed, the solutions do not even conserve atom number. This 
clearly suggests that any explanation of spontaneous emission as being due to 
vacuum fluctuations will not result in an accurate description of the dynam- 
ics. This is a demonstration of the failure of the truncated Wigner method 
for a system which is noticeably simpler than some of those for which it has 
previously been shown to give misleading results [23,24]. 



6 Conclusions 

Using the example of spontaneous emission into a zero temperature reser- 
voir, we have shown how stochastic equations may be developed to model the 



da 
~dt 
da + 



dfi+ 



dt 



dt 

dfi 

~dt 




(34) 
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Fig. 6. Excited state probability for an initially excited one-atom coherent state. 
The solid line is the positive-P average of 9.9 x 10 5 trajectories, the dash-dotted line 
represents 3 x 10 5 trajectories of Eq. 34, and the dotted line is the solution given in 
Eq. 20. 

interaction of bosonic atoms with the electromagnetic field. This approach al- 
lows for the straightforward inclusion of the atomic quantum states, different 
numbers of atoms, the spatial dependence of atomic ensembles and interac- 
tions between the atoms, most of which would be difficult in the usual master 
equation approaches. We have shown that, except in the very special case of 
a one-atom Fock state, the decay is not exponential, but depends on correla- 
tions between the levels as well as bosonic stimulation by the population of the 
lower level. This is manifestation of the many-body atom statistics, because 
the other possible reason, superradiation, is eliminated by our choice of the 
model Hamiltonian. Our approach can be extended to describe the dynamics 
of more complicated processes such as, for example, electromagnetically in- 
duced transparency in degenerate gases. In principle it can also be extended 
to degenerate fermionic atoms using phase-space methods which are under 
development. 
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